Phase diagram of patchy colloids: towards empty liquids 
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We report theoretical and numerical evaluations of the phase diagram for patchy colloidal particles 
of new generation. We show that the reduction of the number of bonded nearest neighbours offers 
the possibility of generating liquid states (i.e. states with temperature T lower than the liquid- 
gas critical temperature) with a vanishing occupied packing fraction (<f>), a case which can not be 
realized with spherically interacting particles. Theoretical results suggest that such reduction is 
accompanied by an increase of the region of stability of the liquid phase in the (T-<p) plane, possibly 
favoring the establishment of homogeneous disordered materials at small (j>, i.e. stable equilibrium 
gels. 



The physico-chemical manipulation of colloidal parti- 
cles is growing at an incredible pace. The large freedom 
in the control of the inter-particle potential has made it 
possible to design colloidal particles which significantly 
extend the possibilities offered by atomic systems [lj. 
An impressive step further is offered by the newly de- 
veloped techniques to assemble (and produce with sig- 
nificant yield) colloidal molecules, particles decorated on 
their surface by a predefined number of attractive sticky 
spots, i.e. particles with specifically designed shapes and 
interaction sites 0- 0. El 0] These new particles, thanks 
to the specificity of the built-in interactions, will be able 
not only to reproduce molecular systems on the nano and 
micro scale, but will also show novel collective behaviors. 
To guide future applications of patchy colloids, to hel 
designing bottom-up strategies in self-assembly H [ 
and to tackle the issue of interplay between dynamic ar- 
rest and crystallisation — a hot-topic related for example 
to the possibility of nucleating a colloidal diamond crys- 
tal structure for photonic applications [9j — it is crucial 
to be able to predict the region in the (T - <f) plane in 
which clustering, phase separation or even gelation is ex- 
pected. 

While design and production of patchy colloids is 
present-day research, unexpectedly theoretical studies of 
the physical properties of these systems have a longer his- 
tory, starting in the eighties in the context of the physics 
of associated liquids [ll El El El El El • 

These stud- 
ies, in the attempt to pin-down the essential features 
of association, modelled molecules as hard-core particles 
with attractive spots on the surface, a realistic descrip- 
tion of the recently created patchy colloidal particles. A 
thermodynamic perturbation theory (TPT) a ppr opriate 
for these models was introduced by Wertheim [l6| to de- 
scribe association under the hypothesis that a sticky site 
on a particle cannot bind simultaneously to two (or more) 
sites on another particle. Such a condition can be natu- 
rally implemented in colloids, due to the relative size of 
the particle as compared to the range of the sticky inter- 
action. These old studies provide a very valuable starting 
point for addressing the issue of the phase diagram of this 
new class of colloids, and in particular of the role of the 



FIG. 1: Schematic representation of the location of the 
square-well interaction sites (centers of the small spheres) on 
the surface of the hard-core particle. Sticks between differ- 
ent interaction sites are drawn only to help visualizing the 
geometry. 



patches number. 

In this Letter we study a system of hard-sphere par- 
ticles with a small number M of identical short-ranged 
square- well attraction sites per particle (sticky spots), 
distributed on the surface with the same geometry as the 
recently produced patchy colloidal particles 4] . We iden- 
tify the number of possible bonds per particle as the key 
parameter controlling the location of the critical point, as 
opposed to the fraction of surface covered by attractive 
patches. We present results of extensive numerical simu- 
lations of this model in the grand-canonical ensemble E3 
to evaluate the location of the critical point of the system 
in the (T - 4>) plane as a function of M. We complement 
the simulation results with the evaluation of the region 
of thermod yna mic instability according to the Wertheim 
theory [la. IL&L Uv^ . Both theory and simulation confirm 
that, on decreasing the number of sticky sites, the critical 
point moves toward smaller <p and T values. We note that 
while adding to hard-spheres a spherically symmetric at- 
traction creates a liquid-gas critical point which shifts 
toward larger <j) on decreasing the range of interaction, 
the opposite trend is presented here when the number 
of interacting sites is decreased. Simulation and theory 
also provide evidence that for binary mixtures of particles 
with two and three sticky spots (where (M) , the average 
M per particles, can be varied continuously down to two 
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by changing the relative concentration of the two species) 
the critical point shifts continuously toward vanishing <j>. 
This makes it possible to realize equilibrium liquid states 
with arbitrary small <j) [empty liquids), a case which can 
not be realized via spherical potentials. 

We focus on a system of hard-sphere particles (of di- 
ameter <7, the unit of length) whose surface is decorated 
by M sites (see Fig. HJ, which we collectively label T. 
The interaction 1/(1,2) between particles 1 and 2 is 
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V(l, 2) = V HS (r 12 ) V w B (r AB ) 



(1) 



AeT BeT 



where the individual sites are denoted by capital letters, 
Vhs is the hard-sphere potential, V^ B (x) is a well in- 
teraction (of depth — uq for x < S, otherwise) and r 12 
and r AB are respectively the vectors joining the particle- 
particle and the site-site centers 20] . Geometric consid- 
erations for a three touching spheres configuration show 

that the choice 8 = 0.5(v5 — 2V3 — 1) ~ 0.119 guaran- 
tees that each site is engaged at most in one bond. With 
this choice of <5, M is also the maximum number of bonds 
per particle. Temperature is measured in units of uq (i.e. 
Boltzmann constant ks = 1). 

To locate the critical point, we perform grand- 
canonical Monte Carlo (GCMC) simulations and his- 
togram re- weighting [2lJ for M = 5, 4 and 3 and for 
binary mixtures of particles with M — 3 (fraction a) and 
M = 2 (fraction 1 — a) at five different compositions, 
down to (M) = 3a + 2(1 ~ a) = 2 A3. We implement MC 
steps composed each by 500 random attempts to rotate 
and translate a random particle and one attempt to insert 
or delete a particle. On decreasing (M), numerical sim- 
ulations become particularly time-consuming, since the 
probability of breaking a bond ~ e 1 ^ becomes progres- 
sively small. To improve statistics, we average over 15- 
20 independent MC realizations. Each of the simulations 
lasts more than 10 6 MC steps. After choosing the box 
size, the T and the chemical potential /i of the particle(s), 
the GCMC simulation evolves the system toward the cor- 
responding equilibrium density. If T and \i correspond to 
the critical point values, the number of particles N and 
the potential energy E of the simulated system show am- 
ple fluctuations between two different values. The linear 
combination x ^ N + sE (where s is named field mixing 
parameter) plays the role of order parameter of the tran- 
sition. At the critical point, its fluctuations are found to 
follow a known universal distribution, i.e. (apart from 
a scaling factor) the same that characterises the fluctua- 
tion of the magnetisation in the Ising model [2l| . Recent 
applications of this method to soft matter can be found 
in Ref. [HEHi. 

Fig. shows the resulting density fluctuations distri- 
bution P{4>) at the estimated critical temperature T c and 
critical chemical potential(s) \i c for several M values |25j. 
The distributions, whose average is the critical packing 
fraction <j> c , shift to the left on decreasing M and become 
more and more asymmetric, signalling the progressive 
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FIG. 2: Density fluctuations distribution P(<j>) in the GCMC 
at the critical point for four of the studied M values. The inset 
shows P(x) for all studied cases, compared with the expected 
distribution (full line) for system at the critical point of the 
Ising universality class [2l|l . 



increasing role of the mixing field. In the inset, the cal- 
culated fluctuations of x, P(x), are compared with the 
expected fluctuations for systems in the Ising universality 
class 21] to provide evidence that (i) the critical point 
has been properly located; (ii) the transition belongs to 
the Ising universality class in all studied cases. The re- 
sulting critical parameters are reported in Table I. Data 
show a clear monotonic trend toward decreasing T c and 
<j) c on decreasing M. 

Differently from the </> c -scale, which is essentially con- 
trolled by M, T c depends on the attractive well width. 
Experimentally, values of uo/ksT comparable to the ones 
reported in Table I can be realized by modifying the phys- 
ical properties (size^polarizability, charge, hydrophobic- 
ity) of the patches U 0, or by furictionalizing the sur- 
face of the particle with specific molecules [2(| H?) • 

One interesting observation steaming from these re- 
sults is that reduction of M makes it possible to shift 
(j) c to values smaller that <p = 0.13, which is the low- 
est <j) c possible for attractive spherical potentials. Indeed 
for spherical square well potentials 0.13 < (j> c < 0.27, 
the two limits being provided by the Van der Waals (in 
which repulsion is modelled by the Carnahan-Starling 
expression) and by the Baxter models, respectively 
with infinite and infinitesimal interaction range. We also 
note that results are consistent with those based on a 
toy model where an ad-hoc constraint was added to limit 
valency and also with previous studies of p articles 
interacting with non-spherical potentials [29ll3C| . 

Visual inspection of the configurations for small (M) 
shows that the system is composed by chains of two- 
coordinated particles providing a link between the 
three-coordinated particles, effectively re- normalizing the 
bonding distance between the M = 3 particles. On 
adding more M = 2 particles, the bonding distance be- 
tween M = 3 particles increases, generating smaller and 
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TABLE I: Values of the relevant parameters at the critical 
point. In the one-component case (M = 3,4,5), fil is the 
critical chemical potential (in units of Mo). In the case of the 
mixture, /xj (jUc) is the critical chemical potential of M = 3 
(M — 2) particles. L indicates the largest box size studied. 
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FIG. 3: Comparison between theoretical and numerical re- 
sults for patchy particles with different number of sticky spots. 
Panel (a) shows the location of the points in the (T - (ft) plane. 
Panels (b) and (c) compare respectively the M dependence 
for T c and <b c . 



smaller <p c . 

To extend the numerical results beyond the point 
where it is currently possible to properly perform GCMC 
(at the lowest (M) each calculation of the 20 studied sam- 
ples requires about 1 month of CPU time on a 3.1 GHz 
processor) and to complement the numerical results, we 
solve the first-order Wertheim TPT El H IH for the 
same model (Eq. ^ . The theory can be applied both to 
one component systems (M = 3, 4, 5) and to binary mix- 
tures ( (M) spans continuously the region from M = 2 
— where no critical point is present — to M = 3) |3lj . 

In TPT, the free energy of the system is written as 
the Helmholtz HS reference free energy Ahs plus a bond 
contribution Ab on d, which derives by a summation over 
certain classes of relevant graphs in the Mayer expan- 
sion |l9j| . The fundamental assumption is that the condi- 
tions of steric incompatibilities are satisfied: (i) no sites 
can be engaged in more than one bond; (ii) no pair of 
molecules can be double bonded. The chosen S guaran- 



tees that the steric incompatibilities are satisfied in the 
present model. In the more transparent (but equivalent) 
formulation of Ref. |32j ] Ab onc i is written as 



(3A, 



bond 



N 



AeT 



Xa 
2 



1 



-M. 



(2) 



Here X A is the fraction of sites A that are not bonded 
The Xas are obtained from the mass-action equation 

1 



X, 



1 + ^2 P x b^ab 

BtT 



(3) 



where p = N/V is the total number density and Aab is 
defined by 



*AB 



= 4?T / 5 , H5(?'l2)(/ J 4s(12)) W i^2 ? 'l2 dr 'l 



(4) 



Here <7ijs(12) is the reference HS fluid pair corre- 
lation function, the Mayer /-function is /as (12) = 
exp{-V£ B {v AB )/k B T) - 1, and (f AB (r 12 )) M rep- 
resents an angular average over all orientations of 
molecules 1 and 2 at fixed relative distance ri 2 . 

The evaluation of requires an expression for 

9hs(t '12) in the range where bonding occurs. We have 
used the linear approximation 12] 



9Hs{r) 



1 f). 

Jl~4 



(l 



2(1 



-(r-1) 



(5) 



which provides the correct Carnahan-Starling |34j value 
at contact. 

To locate the critical point, we calculate the equation 
of state P(V 7 T) = —d(AHs + -^bond)/dVT and search 
for the T and <fi value at which both the first and the 
second volume (V) derivative of the pressure (P) along 
isotherms vanish. Fig. [3] shows a quantitative compar- 
ison of the numerical and theoretical estimates for the 
critical parameters T c and <j> c . Theory predicts quite ac- 
curately T c but slightly underestimates </> c , nevertheless 
clearly confirms the M dependence of the two quanti- 
ties. The overall agreement between Wertheim theory 
and simulations reinforces our confidence in the theoret- 
ical predictions and supports the possibility that on fur- 
ther decreasing (M), a critical point at vanishing </> can 
be generated. 

TPT allows us also to evaluate the locus of points 
where dP/dVr — 0, which provide (at mean field level) 
the spinodal locus. The predicted spinodal lines in the (T 
- <ft) plane for several M values are shown in Fig. 0] On 
decreasing M also the liquid spinodal boundary moves 
to lower <p values, suggesting that the region of stabil- 
ity of the liquid phase is progressively enhanced. It will 
be desirable to investigate the structural and dynami- 
cal properties of such empty liquids by experimental and 
numerical work on patchy colloidal particles. 

We note that our predictions are relevant to a larger 
class of functionalized particles, when particle-particle in- 
teraction is selective and limited in number. Very new 
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FIG. 4: Spinodal curves calculated according to TPT for the 
studied patchy particles for several M and (M) values. 

materials belonging to this class are the recently synthe- 
sised DNA-coated particles ,26]. In this case, M can be 
varied by controlling the number of strands and the at- 
tractive strength can be reversibly tuned by varying the 
length of the strands. Ratios of ua/ksT comparable to 
the ones discussed here can be realized. Again, the phase 



diagram of these new materials has not been experimen- 
tally measured yet and we hope our work will provide a 
guideline. 

For particles interacting with attractive spherical po- 
tentials, phase separation always destabilises the forma- 
tion of a homogeneous arrested system at low T. Instead, 
it is foreseeable that, with small (M) patchy particles, 
disordered states in which particles are interconnected in 
a persistent gel network can be reached at low T with- 
out encountering phase separation. Indeed at such low 
T, the bond-lifetime will become comparable to the ex- 
perimental observation time. Under these conditions, a 
dynamic arrest phenomenon at small <fi could take place. 
It would be possible to approach dynamic arrest continu- 
ously from equilibrium and to generate a state of matter 
as close as possible to an ideal gel |3fij . 

The study of the structural and dynamic properties of 
these low (M) equilibrium systems will hopefully help 
developing a unified picture of other interesting network 
formation phenomena taking place at low <f> [3a . Wh . l38l 
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